AAS 04-287 


REPRESENTATIONS OF INVARIANT MANIFOLDS FOR 
APPLICATIONS IN THREE-BODY SYSTEMS 

K. Howell,* M. Beckman/ C. Patterson,** and D. Folta 5 

The Lunar Li and L 2 libration points have been proposed as gateways gr antin g inexpensive 
access to interplanetary space. To date, only individual solutions to the transfer between three- 
body systems have been found. The methodology to solve the problem for arbitrary three-body 
systems and entire families of orbits is currently being studied. This paper presents an initial 
approach to solve die general problem for single and multiple inpulse transfers. Two different 
methods of representing and storing the invariant manifold data are presented. Some particular 
solutions are presented for two types of transfer problems, though the emphasis is on 
developing the methodology for solving the general problem. 

INTRODUCTION 

With the increasing interest in missions involving Sun-Earth and Earth-Moon libration points, 
it is necessary to further develop numerical, and possibly semi-analytical, tools to assist in trajectory 
design in multi-body regimes, including libration point orbits. Halo orbits are well-known examples of 
periodic orbits in such regions of space. Thus far, Earth-to-halo transfers, as well as halo-to-Earth arcs, 
have been computed using a number of different numerical procedures including exploitation of the 
invariant manifolds associated with a particular periodic halo orbit (or quasi -periodic Lissajous 
trajectory). More recently, transfers between different three-body systems are a new focus for potential 
mission scenarios. Developing transfers between Earth and a halo orbit or between different three-body 
systems involves numerically integrating trajectories from different initial conditions near a desired 
halo orbit manifold until one is identified feat is most suitable for the application of interest. Of 
course, the stable/unstable invariant manifolds that correspond to a single halo orbit are infinite in 
number, but reside on the surface of a tube. 1 Nevertheless, continuous computation of individual 
manifolds using numerical integration is not efficient or even practical for some applications. In the 
analysis of trajectories to/from a halo orbit, for example, the size of the most useful periodic orbit may 
be unknown and its amplitude may serve as a design parameter. The design space then includes not 
just a tube corresponding to the invariants manifolds of one halo orbit; rather it becomes a volume 
consisting of many tubes. For the problem of system-to-system transfers, the goal is the intersection of 
two manifold tubes - one from each system. A maneuver at an intersection point will shift the vehicle 
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from one tube to the other. In the Circular Restricted Three Body Problem (CR3BP), the flow in this 
region of space can also be visualized by noting that the tube is a separatrix, i.e., it bounds different 
regions of the flow. Thus, transit trajectories might be sought that pass inside the tubes and shift from 
one tube to the next. In this case, the spatial intersection of the tubes is still the key information to 
complete the computations. Computing those intersections is facilitated by efficient but accurate 
approximations of the tubes. So, the objective of this work is twofold: (i) 

representations/approximations of the invariant manifolds associated with periodic libration point 
orbits that are efficient and accurate in an automated process; (ii) demonstration that the approximation 
can be used to generate solutions to representative problems. 

It is generally impossible to determine an analytical expression for stable and unstable 
manifolds so accurate numerical computations are required. This method usually proceeds by using 
one set of initial conditions and numerically integrating a single trajectory. Increasingly, the required 
trajectory length is long and accuracy suffers if the step size is small enough to accommodate 
potentially sharp changes along the path. Hobson offers some consideration of this issue and a method 
for approximating the manifold that includes an estimate of the errors. 2 Hobson’s method is an 
improved approach to globalize the manifold in terms of step size, but it still employs the standard first 
order approximation to initiate the process. A single trajectory may not, however, reflect the actual 
complex dynamical behavior of the system. In fact, the problems of interest here include a complete set 
of manifolds associated with a periodic orbit as part of the potential solution space. In configuration 
space, the stable or unstable manifolds are located along the surface of a tube corresponding to the 
periodic halo orbits. Either the numerical computation of some minimum number of orbits is required 
to represent the tubes or individual trajectories must be computed or recomputed as needed. If such 
data must be available for analysis, it may require a large amount of storage and retrieval capability. 
The specific requirements depend upon die application, of course. Relatively recently, a number of 
mathematicians have investigated approximation techniques. Guder et al. 3 introduce cell mapping for 
the prediction of long term behavior and global analysis of nonlinear dynamical systems. Dellnitz and 
Hohmann 4 use these cell mapping techniques combined with a subdivision procedure to approximate 
unstable manifolds and global attractors. However, they are not interested in determining the foil 
global behavior; rather a numerical method that allows the approximation of the global attractor up to a 
specified accuracy. A box in IR" is specified in which the dynamical behavior is to be analyzed. The 
box is subdivided and boxes discarded that do not contain part of the relative global behavior. A 
continuation method is further developed by Dellnitz and Hohmann to approximate the types of 
invariant manifolds of interest here. This method is also called an “outer approximation” of the 
unstable manifold. It is very useful to obtain a global view, however, it does not offer a 
parameterization of the manifold via arc length and small details may cause difficulty. Another 
technique by Krauskopf and Osinga 5,6 is designed to specifically produce a list of points on the 
manifold by adding new points at prescribed distances from the last point and essentially “growing” 
the manifold. It is difficult to search this type of solution space, however, for the state information that 
is desired for the applications presented in this paper. 

Approximating the manifold tubes for application to mission design is the first step in this 
current investigation. The process must be straightforward and offer accurate representations of the 
position and velocity states. It might also be used to represent the volume that includes multiple 
manifold surfaces to aid in determining a specific halo orbit for a given scenario. This process is 
applicable to error analysis, recovery and design for contingencies. Recent studies suggest that these 
types of analyses would greatly benefit from another parameter such as halo orbit size; if the manifold 
is represented in an easily accessible form, it can be accomplished in a automated manner. 7 ’ 8 
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Application of these approximations to the system-to-system problem of recent interest is also 
used to demonstrate the approach. As noted by Farquhar et al., 9 the libration points may become a 
primary hub for future human activities in the Earth’s neighborhood. The Sun-Earth L 2 point is 
expected to be where a number of large astronomical observatories will be located. These L 2 telescopes 
may require human servicing and repair as the missions and hardware become more complex and 
expensive. Closer to Earth, the Earth-Moon Li libration point may serve as the staging node for the 
missions to the Sun-Earth L 2 point as well as the Moon, Mars, and the rest of the solar system. Folta et 
al. 10 produce various trajectories to support the human servicing role and offer other examples as well. 
But, improved and automated schemes for computation of such designs will be necessary for efficiency 
and larger studies. Lo and Ross" also support the Earth-Moon Li point as a “portal” to move beyond 
the Earth’s neighborhood. The manifold tubes are introduced by Lo and Ross as the basis for the 
design strategy to produce the trajectories to move between these systems. Koon et al. 12 describe the 
concept in more detail for the Earth-Moon to Sun-Earth problem. It has also been examined in the 
Jovian system. 13 In any case, intersections may ultimately be sought between many tubes from many 
different halo orbits in each system; the complexity forces a new look at the computations. Of course, 
this problem creates a very complicated solution space and some representation of this set of solutions 
is sought to facilitate an automated design procedure. This study examines both splines and a cell 
structure to represent the manifolds. 

MANIFOLD APPROXIMATIONS: SPLINES 

In order to consider the intersection of a single trajectory (or a group of trajectories) with a 
surface that defines an invariant manifold corresponding to some halo orbit, the manifold must be 
generated and stored in an accessible format. One method investigated is the use of splines. The 
existing algorithms within the MATLAB™ Spline Toolbox™ are incorporated. 

The numerical manifold data set that corresponds to one periodic halo orbit in some arbitrary 
system is first generated using a variable step propagator and the differential equations that model the 
circular restricted three-body problem (CR3BP). Multiple individual trajectories along one tube are 
generated to produce a good approximation to the entire manifold surface. For this analysis, N points 
along the halo orbit are identified and then globalized to generate N manifolds identified with one 
tubular surface corresponding to a specified periodic halo orbit. The initial state vectors are propagated 
and state vectors are identified along each trajectoiy, each at the same time, that result in Nx 6 vector 
elements that together represent a manifold surface. This data set is placed in a grid for later 
computations. The grid is comprised of A columns and M rows (M number of integration points, each 
at the same time). Spline knots are placed along each trajectoiy at the location of each isolated state 
vector, measured in time from the initial state at the halo orbit. The spline function is generated from 
the knot sequence, ^ , and the polynomial coefficients, c jt , 

Pj(*)= Eti c ji J = \:n (1) 

A bivariate manifold spline is then generated for a manifold surface corresponding to a single halo 
orbit. The bivariate spline function yields manifold position data for a given trajectory, rj, or 
equivalently the phase angle or time of the originating halo orbit (a total of A trajectories), and the time 
from the start of the manifold, r (a total of M points), 
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This process enables calculation of position components given two parameters: number of the 
trajectory (or trajectory tag number) and time from the periodic orbit. The result appears in Figure 1 . 



Sun-Earth L, (SEL1 ) Manifold Earth-Moon L 2 (EML2) Manifold 

Figure 1 Position Splines 


MANIFOLDS APPROXIMATIONS: CELLS 

Also investigated is a representation of the manifolds in terms of a relatively 
straightforward cell structure such that each cell includes a polynomial model of one or more 
manifolds. Thus, the potential set of solutions may incorporate manifolds from various halo 
orbits. The cell structure then includes volume elements that can be searched for intersections 
with other single trajectories or intersections with other cells. 

In the Sun-Earth system, consider a stable manifold tube that is computed numerically. 
With this initial surface, the manifold tubes that correspond to several additional, nearby halo 
orbits fill a volume in space as the individual manifold surfaces wrap around each other. To 
model any particular tube, exploit the fact that over small regions the stable manifold of a halo 
orbit in the circular restricted three-body problem is nearly flat and can be represented in position 
and velocity by low order functions. A volume of space through which several manifold surfaces 
pass and that is sufficiently small to be approximated in this maimer is denoted here as a cell. 
(See Figure 2.) 


Figure 2 Cells are Volumes in Space Defined to Contain 
Small Regions of Manifold Trajectories 

Within a cell, position data obtained from a single stable manifold tube are approximated using 
the fit function from Mathematica® to produce an analytical function of a surface. This process is 
repeated with other nearby manifolds until surface approximation functions for several manifolds 
passing through the cell are available. With this same position data, velocity data are also used to 
generate a single fit of each velocity component within the cell as a function of position . Thus, if 
N stable halo manifolds pass through a single cell, N fits to approximate the manifold surfaces are 
produced as well as Nx3 fits to approximate velocity components as a function of position. 

Fitting and Cell Shapes 

The velocity components are fit as third order polynomial functions of position in the 
standard rotating coordinate frame associated with the circular restricted three-body problem 
(CR3BP), that is, 

v w ( x > y*- z ) - 1 a o +• °v x + a i x 2 +i a r * 3 + a 4 y + a s x y +■ q 6 * 2 y + 

a 7 y 2 +a i xy 2 +a 9 y 2 +a ]0 z + a u xz + a n x 2 z + (3) 

a n yz + a u xyz +a, 5 y 2 z + a 16 z 2 + a 17 xz 2 + a xi yz 2 + a l9 z 3 

For the surface (position) fits, which often globally exhibit significant curvature, the position data 
are transformed to a spherical coordinate set and the surface is defined by a fit of radial data as a 
function of angular coordinates as follows, 

r(0, <t>)=a 0 + a, cos(^) + a 2 cos 2 {<f>) + a 3 cos(0) +a 4 cos{<f>) cos(0) + 

a s cos 2 ( 0 ) + a 6 sin(^) + a 7 cos(^) sin(^) + a % sin(^) cos(0) + 

2 ( 4 ) 

a 9 sin {<f>) +a ]0 sin(0) + a, , cos(^) sin(0) + a 12 cos(f?) sin(0) + 

a, 3 sin(0) sin(0) + a 14 sin 2 (0) 


where 9 and (p are angles that orient the radial direction in 3D. The origin for the spherical set of 
coordinates is selected as the average center of curvature computed from every data point within 
the cell. This transformation is accomplished by defining unit vectors in the average tangential, 
normal, and binormal directions from the mean center of the data and placing the cell origin along 
the normal direction at a distance equal to the average radius of curvature. With a small data set, 
this choice encourages a cell shape that is indicative of the local shape of the manifold itself and a 
smooth surface fit with little variation in r. In some areas of the tube, however, manifold 
trajectories tend to lie in the same plane as the path representing the centers of curvature for the 
data; for a cell to contain several such trajectories this choice will result in a large variation in r. 
In this case, the origin is shifted to a location along the binormal direction. With the data points 
transformed to this frame, the cell boundaries are defined as the minimum and maximum values 
of r, 6 , and <f > . The fits are reasonably accurate everywhere except very near the Earth where 
curvature of the actual manifold is extremely high. Along the manifolds, a high curvature region 
requires the use of smaller cells within which the manifolds can be approximated as nearly flat. 
Also, the volume of space contained in a cell should not extend too far beyond the volume 
defined by the manifold. Consequently, cell sizes and shapes are tailored to the local shape of the 
manifold volume. Then, position and velocity functions for the manifold data within the cells can 
be determined with standard fitting algorithms. This process has been automated for examples of 
the type examined here. In Figure 3(a), the maximum and median position errors along the 
surfaces are compared with actual manifold tube data within the cells. Median errors are 
typically less than 500 km, and maximum errors typically less than 5000 km, the most notable 
exceptions being very near the Earth where the error actually extends beyond the range of the 
display. In Figure 3(b), the median and maximum percentage errors in the velocity functions are 
compared to actual velocity data. Note that the velocity errors are typically very small. 
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Figure 3 SEL1 Manifold Position and Velocity Fits :Halos Range in Size from 
Az « 800,000 km to Az « 880,000 km 
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Cell Creation 


Approximation of an entire volume of manifolds involves slicing that volume into many 
such cells. The creation of the cells begins with a single tube of trajectories, each propagated 
backwards to its first periapsis point. The cells cover the entire area of the tube and are kept 
relatively small so that the analytical approximations remain accurate. Also, they are designed to 
overlap at their edges so that the fits within one cell mesh smoothly with the fits of a neighboring 
cell. Subsequently, the fits within the cells collectively represent a more continuous 
approximation of the manifold. 

To create the cells, a manifold tube is sliced into two sections. The first section is formed 
by truncating the tube close to the periodic orbit. Using the libration point as an origin, each 
revolution along the path as it departs the orbit is sliced typically into 64 sectors and the data 
within each sector are used to define a new cell and frame as described above. This process is 
then repeated a predetermined number of times. The second section is comprised of the rest of the 
tube. The method of cell creation is then modified to more closely follow the evolution of the 
manifold trajectories. A limited number of neighboring trajectories on the manifold, denoted here 
as a “ribbon”, is propagated backward from the last section to periapsis and sliced evenly in 
distance along its length. Each slice of data taken from the ribbon is used to define a new cell. 
The process is repeated using the last trajectory of the previous ribbon as the first trajectory of a 
new ribbon until the entire the tube is represented. Note that "before this step is complete, some 
cells may require resizing to maintain accuracy. A maximum volume is selected and any cell 
larger than the maximum is sliced in half along its largest side. When placing data in cells it is 
also important that a minimum amount of data from a minimum number of trajectories always be 
present so that surface fits can be accurately made. To accomplish this, manifold tubes have 
actually been propagated twice, the second time containing three times the number of trajectories 
than the first, and with many more points. The manifolds with the sparse data are used for most 
of the cells created, but when the minimum point and trajectory requirements are not met within a 
cell, the densely populated arc is incorporated. 

Adding Manifold Data 

After initial cell creation, the cells contain the data corresponding to merely one manifold 
tube. Data from additional manifold tubes can be added with certain caveats. The cell shapes at 
this point are tailored to the original tube shape and size, so the trajectories of any other manifold 
tube will not be fully contained within their volumes. Therefore, to contain as much data as 
possible from a nearby tube, the cells may be required to increase in size. However, the fitting 
process is premised on the idea that manifold data can be approximated only over small regions, 
so cell growth is limited to encompass only the new tube data. A single cell originally contains all 
manifold data from one tube within a bounded volume. The data of a second tube is transformed 
to the original cell’s spherical coordinate frame, and the angular bounds of the cell are increased a 
small amount. Any part of the new tube’s data that fits within the new angular bounds of the cell 
is considered as data for a second surface within the cell. The bounds of the cell are then reset to 
reflect the increase in the volume of data stored. In Figure 4 is an example of three SEL1 
manifold tubes which have been placed in the same set of cells. They range in size from 
Az=~800,000 km to Az =-880,000 km. It is apparent that even with three tubes placed in the 
cells, there can be vacant areas. In these areas, the velocity fits will be less accurate. 
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This problem can be solved by adding more 
manifold tubes without allowing cell growth, but 
maintaining the minimum number of points and 
trajectories. If, in a single cell, the minima are 
not met then the data for that surface does not 
get added to the cell. So some cells may contain 
more manifold tubes from a wider range of halos 
than other cells. It is also possible to 
approximate the wider range of halos by 
stacking cells, essentially implementing the cell 
creating process for the manifolds of one set of 
halos, then creating entirely new cells for a 
larger or smaller set and repeating. Figure 4 Five Tubes through 

a Set of Cells 



Example: Using Cells to Find Halo Transfers in the CR3BP 

With approximations defined for stable Sun/Earth Li (SEL1) halo manifolds within the 
cells, the process of finding transfers to halo orbits is simplified. For any spacecraft trajectory 
that intersects a cell at a point r mt with a velocity v jm , an initial guess for a AV that might be 

necessary to complete a transfer is simply AV = V^,)- where V^) is the approximation 

for a manifold velocity at the intersection point. It is assumed that the continuous distribution of 
periodic orbits in a family possesses manifolds that exist on a continuum of surfaces passing 
through the cells, and that velocity on those surfaces varies continuously throughout the volume. 
Analytical approximations for only a few of the surfaces of position are planned to be available, 
but a single velocity approximation is continuous over the entire cell volume. So, even if Tin, lies 
between or outside the known manifold surfaces within the cells, it is still identified with some 
member of the manifold continuum, thus, the velocity fit can be applied and used as an initial 
guess for the velocity of a manifold trajectory passing through that point. Moreover, it is not 
necessary that a particular halo orbit be specified in advance. Indeed, if the best intersection point 
within the cell is to be selected from among all possible points in the volume, and that best point 
is not known a priori, then any halo orbit in or near the range of those used to initially generate 
the manifolds and create the cells may serve as the ideal target halo for the transfer. So rather 
than selecting a halo and attempting to locate intersections with its stable manifold, an easier 
approach may be the determination of the intersections with a volume of stable manifolds that are 
associated with a set of periodic orbits and use the lowest-cost intersection point to define the 
target orbit. 

The Use of Surface Fits in Determining a Destination Halo 

An intersection point r^, inside a cell has spherical coordinates This point may 

be inside, outside, or somewhere between the surfaces defined within the cell that correspond to 
known manifold tubes as is apparent in Figure 5. To discover its location relative to the existing 
surface fits, calculate rfO,#) where i ranges from 1 to n, the number of surfaces within the cell. 

Each surface, and thus each value of is associated with a specific halo. These halo 
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To halo #1 



orbits can be characterized by the initial 
conditions used to propagate them 
numerically in the 3-body problem. For 

example, (x' o ,0,z' o ,0, j^O) are the states 

at the point of the maximum z amplitude 
and are used to parameterize a single halo 
orbit. Each initial condition is fit as a 
quadratic function of ^(9,0), directly 
associating the halos to the sizes of the 
manifold tubes at {0,<j)) within the cell. 
This procedure yields functions that can 
be solved for a set of initial conditions 
(and thus the size) of a destination halo 
orbit as a function of the location r within 
the cell, that is, 


x Q =a x0 + a x s + a x2 r 2 

(5a) 

Zo=a Z Q+a zX r + a z2 r 2 

(5b) 

To = a 9Q + a n r + a n r 2 

(5c) 


Figure 5 Locating the Intersection Point within a Cell 


The coefficients can be determined by solving the following three matrix equations. 
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If the number of surfaces n > 3, this is accomplished with a least squares fit. If n < 3, then a less 
accurate linear fit must suffice. However, an effort has been made to always include at least three 
manifold surfaces in eveiy cell. 

Approaching the Destination Halo 

Approximations from fit functions are generally most successful when applied in a region 
near or within the bounds of the region where the data are available. This region is defined by the 
surfaces within the cell. Thus, the surfaces nearest and furthest from the cell origin define the 
region. The distance between these two surfaces at any point (9,<j>) in the cell defines a width, 
w(6,<f>), and any trajectory intersection point that is more than a distance w from the surfaces is 
ignored as a potential transfer point for the given cell (but, it would re-emerge in another cell). 
These low order velocity fits within the cells are generally sufficient to produce a trajectory arc 
from a point within the cell to the vicinity of the libration point orbit, but it is likely that it will not 
directly insert into a periodic orbit. The sensitivity of the problem demands that some amount of 
targeting is used to complete the transfer. However, with the approximation, it is an additional 
automated step that results in successful transfers to halo orbits in the circular restricted problem. 

Assume that an initial trajectory originates from somewhere in or beyond the current 
system. There is no restriction on the source of the incoming leg. Assuming that this path will 
intersect one or more cells, they are identified and all intersecting points rim within the cell are 
tagged. For each intersection point, the velocity discontinuity is calculated, that is, 

AV cgS = V (r^, ) - Vjjj. seeking the r int that minimizes AV,*,,. The calculated AV«n is applied. 

The point of closest approach to the halo is determined. Since AV ceU is an approximation, the 
transfer arc will not approach the periodic orbit asymptotically. A halo orbit insertion (HOI) point 
is selected and targeted from rj„, and a HOI AV is included at the HOI point. This completes the 

transfer with a total cost C = |MV ce „ + AV (ar?e(er | + (AV^ | This cost is a function of the HOI 

point. For convenience, the MATLAJB™ Optimization Toolbox™ function fminbnd is utilized to 
find a minimum of this cost as a function of the HOI location. 


Numerical values for a test case appear in Table 1 and Figure 6. Three halos at Az = 
-800,000 km, Az = -840,000 km and Az=~880,000 km are used to generate manifold cell 
volumes, and a fourth stable manifold to a halo of Az = -945,000 km, one that is not used in any 
of the fitting, is used to test the transfer procedure. A point along this actual stable manifold is 
perturbed with a AV of 20/ms at a point within a cell. At this precise 


Table 1 

Targeting a Halo Orbit 


Step in Process 


Exact Manifold 
Intersection 
Approximation Result 
at Same Intersection Point 


Destination 

Halo 

Az = 945,48 lkm 
Az = 945,176 km 


Total AV at AVhdi 
intersection 
20 m/s 

19.69 m/s 1.43 m/s 
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intersection point, the approximations within the cells are employed to estimate the velocity at 
this point, approximate the appropriate AV and determine the specific halo that would minimi?^ 
the cost. The actual manifold should be the best solution and die results of the approximations 
appear in the table. A good estimate of the destination halo is achieved and the transfer is 
accomplished with a cost veiy similar to the 20m/s perturbation. 


X ID 5 

HOI 



Figure 6 Transfer Determined Using Cells 


SYSTEM TO SYSTEM TRANSFERS 

A more challenging problem is the system to system transfer. In particular, for 
demonstration here, the transfer from the Earth-Moon system to the Sun-Earth system is 
considered. As noted previously, these transfers are the focus of a number of researchers. 14 In this 
work, an efficient method of computation is sought using approximations to the manifolds. The 
ultimate objective is to use the approximation to produce the transfers in full, 3-D ephemeris 
models that may include other perturbations as well. This problem is approached using both 
splines and cell structures. Ultimately both may be useful in designing these types of trajectories. 

Using Splines to Determine the Intersections of Earth-Moon and Sun-Earth Tubes 

Using the manifold position spline functions, intersections can be computed between 
invariant manifolds of different three-body systems. Assuming that the goal is to shift directly 
from the surface of one manifold to the other, these intersections result in one-impulse transfers 
between the three-body systems. Such could be the case if it was desired to depart an Earth-Moon 
L 2 orbit and reach a Sun-Earth Lj halo orbit, for example. Figure 7 includes an example of a large 
Sun-Earth Li stable manifold and various Earth-Moon L 2 unstable manifolds at different phases 
of the Moon. From this figure, it is apparent that the tubes intersect. Using a 
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SEL1 manifold 



X (km) Y (km) 


Figure 7 Intersections Between Tubes 


Moon phase of 210 deg, an intersection between the manifolds is assured at greater than one 
point. The Matlab Optimization Toolbox is employed to determine the intersections. The Matlab 
function fmincon yields a minimum of a constrained nonlinear multivariate function. Medium- 
scale optimization is used and the gradient is determined through finite-differencing. The SQP 
method closely mimics Newton's method for constrained optimization just as is done for 
unconstrained optimization. At each major iteration, an approximation is made of the Hessian of 
the Lagrangian function using a quasi-Newton updating method. This Hessian is then used to 
generate a quadratic programming subproblem whose solution is used to form a search direction 
for a line search procedure. The solution yields a position on each manifold, the two of which 
differ by only 43 km for this example. The first successful use of fmincon gives us one point of 
intersection. It is necessary to find additional points of intersection where an entire set should be 
closed ring (assuming complete intersection of entire tube). Such a ring can be computed by 
adding additional constraints to fmincon. A nonlinear equality constraint forces fmincon to locate 
a second point of intersection between the two manifolds but at a distance h from the first point. 
The distance h is currently set at 1000 km as an initial guess. The ring is successfully computed 
and appears in Figure 8. Once the intersection is determined in configuration space, the velocity 
differences are also available, and the cost can be computed. 
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Figure 8 Ring of Intersection 


Using Cells to Approximate a Transfer and Shift into the Ephemeris Model 

The same Earth-Moon system to Sun-Earth system transfer was also attempted using the 
cell structure. In this case, the problem is solved in two steps. First, an approximate solution is 
obtained in the circular restricted problem using the cells. Then, the solution is transferred to a 
more complete model. Using the software package Generator, developed at Purdue University, 
this approximation is used to produce an end-to-end trajectory from an Earth-Moon libration 
point orbit to a Sun-Earth libration point orbit. It is noted that other researchers are also seeking 
transfers of this type in the full ephemeris model. 10 ' 12 

In the circular restricted Earth-Moon system, 
unstable manifold trajectories for a Li halo of 
amplitude Az = -39,000 km are propagated. These 
manifolds form a tube, of course. Any of the 
trajectories can be transformed to Sun-Earth 
coordinates and intersections are determined with 
cells that correspond to the SEL2 stable manifold. A 
set of cells is generated for SEL2 stable manifolds 
to halos ranging from Az = -170,000 km to Az = 
-220,000 km as seen in Figure 9. These cells 
contain position and velocity fits for the data of 

only three tubes. Initially, the Moon is 45° out-of- 

Figure 9 EML1 Paths Intersect Cells 
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phase with the Sun-Earth system. It is also noted that the EM frame is rotated -5° about the Sun- 
Earth x-axis to more accurately represent the inclination of the Moon’s orbit. After the 
transformation to the Earth-Sun system, the Earth-Moon manifold trajectories are searched 
locate all intersections with the SEL2 cells, and the transfer approximation process proceeds. 
Figure 9, all EM manifold trajectories are overlaid with the SEL2 cell outlines. The estimated 
total cost for the “best” transfer is 24.518 m/s, as determined in the circular restricted problem. 
The transfer appears in Figure 10 in both the circular restricted Earth-Moon view (Figure 10(a)) 
and the circular restricted Sun-Earth view (Figure 10(b)). In the Earth-Moon view it is observed 
that the manifold that departs L] does encircle L 2 as it departs the system. This characteristic is 
actually familiar. 



(a) (b) 


Figure 10 Transfer from the Earth-Moon L, Halo Orbit to a Sun-Earth L 2 Halo Orbit 


Table 2 

Solution for Transfer to Halo 

Destination Total AV at AVh™ 

Halo intersection 

Approximation Result Az= 123,265 km 21.146 m/s 3.372 m/s 

at Best Intersection Point 
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ET S 





Recall that the AV here is calculated in the Sun-Earth frame with no regard for the dynamic 
effects of the Moon. This solution is transferred into a model with a complete ephemeris 
formulation. The Generator 3.0.2 software, developed at Purdue, incorporates an RK 8/9 
integrator; the gravity of the Sun, Earth, and Moon is included for this example as well as the 
modeling for their ephemeris locations, but no solar radiation pressure force is added. The first 
arc of the transfer, as viewed in the Earth-Moon frame, appears the same as that in Figure 10(a). 
The Sun-Earth view is seen in Figure 11 and some differences relative to the solution in the 
circular problem now appear. The cost at the intersection point is now 15.27 m/s; there is no 
longer any libration orbit insertion AV . 


x 10* 



x lQ r ' 


x ID 3 



x 10* 


Figure 11 Transfer in the Sun-Earth System: Ephemeris Model 


Finally, patch points are introduced into Generator (as marked in Figure 12) and all legs are 
simultaneously corrected to yield a zero cost solution. Of course, the final libration point orbits 
are not periodic but are Lissajous quasi-periodic trajectories. The final size of the halo orbits in 
the full model in the Earth-Moon system corresponds to Az = 39,600 km; in the Sun-Earth system 
the out-of-plane amplitude is Az = 163,000 km. In the Sun-Earth view, the final ephemeris zero 
cost solution cannot be distinguished from the plot in Figure 1 1. It is interesting to compare the 
Earth-Moon departure path in the approximation in the circular problem to the final ephemeris 
corresponding to the no cost solution. This comparison appears in Figure 13. All of the general 
characteristics are maintained. As this example is reviewed, observe that every step of the process 
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yields a libration point orbit of different size. This flexibility appears to aid the search for a 
solution. 


x 10* 



Figure 12 Patch Points in the Cost Reduction Algorithm 


x io 5 



x 10 s 


Figure 13 Comparison of the Circular Approximation to the Departure trajectory from the 
Earth-Moon System with the Final Ephemeris Solution 
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SUMMARY 


The objective of this research activity is to investigate approximations for stable and 
unstable manifolds to use these objects in mission design applications. Both splines and cell 
structures are investigated and prove successful in computations. Further developments are 
planned that include a wider range of applications. 
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